\chapter{寻找函数极值}

\subsection{一维函数最优化}
有时候一个物理量f是一个参数x的函数$f(x)$，而我们需要找到最优参数，使得$f$取得最小值（或者最大值）。
这是很常见的数值问题，比如，量子力学中的变分法，就是一个求极小值的问题。
最简单的方法是，我们在参数的取值范围内取一系列离散的点$x_i$，分别计算其函数值$f(x_i)$，然后取其中的最小值，作为函数在这个区间内的最小值，如图\ref{fig:optimization_demon}所示。
\begin{figure}[h]
\centering
\includegraphics[width=0.7\linewidth]{optimization_demon}
\caption{示意图：寻找函数极值，理论上的极值为x=-0.5。}
\label{fig:optimization_demon}
\end{figure}

如果希望做得更加精确一点，就可以用算法来逼近极小值。
\begin{shaded*}
\begin{itemize}
\item 第一步是找到一个区间，确保这个区间内有极小值，比如图\ref{fig:optimization_demon}中区间$[-1,1]$。
\item 第二步是不断缩小这个区间，逼近极小值。
\end{itemize}
\end{shaded*}

第一步比较容易，可以如下实现
\begin{shaded*}
\begin{itemize}
\item [1] 给定一个初始区间$[a,b]$, 比较$f(a), f(b)$两个值，如果$f(a)<f(b)$，调换$a,b$的值。寻找极小值的方向为$ a \rightarrow b \rightarrow $。
\item [2] 定义$c = b + w*(b-a)$, 其中$w$为一个常系数，代表了每次前进的相对步长大小，有人取作 1.618。
\item [3] 如果$ f(c) < f(b) $，说明函数值还在减小，则令$a=b, b=c$，回到步骤2。
\item [4] 如果$ f(c) > f(b) $，说明$ [a,c] $区域内有极小值。
\end{itemize}
\end{shaded*}

比如，在图\ref{fig:optimization_demon}示例中，若初始区间为$[a,b]=[0,0.1]$，则会在第1步得到$[a,b]=[0.1,0]$，在第2步得到$c=-0.1618$，反复执行几次2,3步骤以后，在第4步得到一个区间，确认该区间内有极小值。

第二步是缩小“嫌疑空间”，逼近极小值。
如果初始区间为 $[a,c]$，其中有一点 $b$，满足 $f(a)>f(b), f(c)>f(b) $，可以如下实现完成第二步

\begin{shaded*}
\begin{itemize}
\item [1] 取 $[a,c]$ 区间内另一点 $d$，一般取 $[a,c]$ 的某个比例分割点，比如黄金分割点$d=a+0.618(c-a)$，求得 $f(d)$ 的值。
\item [2] 如果$d>b$：①如果 $f(d)<f(b)$，则 $a=b, b=d$；② $f(d)>f(b)$，则 $c=d$。
\item [3] 如果$d<b$：①如果 $f(d)<f(b)$，则 $c=b, b=d$；② $f(d)>f(b)$，则 $a=d$。
\item [4] 如果$|a-c|$小于某个指定的精度，就将$(a+c)/2$当作极小值位置。
\end{itemize}
\end{shaded*}

实现上述功能的代码如下
\begin{lstlisting}
#include<iostream>
using namespace std;

#include<cmath>

#define max_iteration 1E4

/*
 * findmin: 给定函数func, 初始区间[a,b]，在精度要求precision下，寻找一个极小值点，并返回该点坐标
 */
double findmin(double (*func)(double), double a, double b, double precision){

        double temp,c,d;
        int i;

        if( func(a)<func(b) ){//寻找极小值的方向为：a->b->
                temp=a;
                a=b;
                b=temp;
        }

        for(i=0;i<max_iteration;i++){//寻找嫌疑区间
                c= b + 1.618*(b-a);
                if( func(c) > func(b) )break;
                a=b;
                b=c;
        }
        if(i==max_iteration){//10000次迭代以后，没有找到嫌疑区间，报错
                cout<<"after "<<max_iteration<<" iterations, didn't find a zone for fmin"<<endl;
                exit(1);
        }

        for(i=0;i<max_iteration;i++){
                if( fabs(a-c)<precision )break;
                d=0.38197*a + 0.61803*c;
                if( (d-a)*(d-b)<0){//如果d在a,b之间
                        if( func(d) < func(b) ){
                                c=b;
                                b=d;
                        }
                        else{
                                a=d;
                        }
                }
                else if( (d-a)*(d-b)>0 ){//如果d在b,c区间
                        if( func(d) < func(b) ){
                                a=b;
                                b=d;
                        }
                        else{
                                c=d;
                        }
                }
                //下面四行帮忙查错：观察每次迭代的结果
                //cout<<"-----------------------------------"<<endl;
                //cout<<"a="<<a<<", f(a)="<<func(a)<<endl;
                //cout<<"b="<<a<<", f(b)="<<func(b)<<endl;
                //cout<<"c="<<a<<", f(c)="<<func(c)<<endl;
        }
        return (a+c)/2;
}

double f(double x){
        return exp(x);
}

int main(){
        double x = findmin(f, 0, 0.1, 1E-6);
        cout<<" fmin is at x="<<x<<endl;
        return 0;
}
\end{lstlisting}

\subsection{共轭梯度法}
\subsubsection{引言}
之前我们讲过最小二乘法，可以优化形如 $f(\vec{x}) = \sum\limits^n_{i=1} a_i g_i(\vec{x})$ 的公式，使之尽可能地接近已有的实验数据点。
%当时我们定义了 $\Delta (\vec{a}) = \sum\limits^N_{j=1} ( f(\vec{x}_j) - y_j)^2$，寻找 $\Delta(\vec{a})$ 的极小值点，来得到最优化参数 $\vec{a} = (a_1, a_2, \cdots, a_n)^\mathrm{T}$。
共轭梯度法是一种更一般的优化方法，可用于寻找多元函数 $f(\vec{x})$ 的极小值点。

在计算中，首先需要猜一个初始位置 $\vec{x}_0$，在这个方向上沿 $\vec{p}_0$ 方向寻找极小值点（即一维最优化）。
这个极小值点是 $\vec{x}_1$，沿一个新的方向 $\vec{p}_1$。
如此进行下去，在 $\vec{x}_i$ 点转弯，沿一个新的方向 $\vec{p}_i$ 寻找极小值。
沿着这样的路径，最终找到极小值点。

也许你会直觉地认为，在每个 $\vec{x}_i$ 点出发，应该沿梯度的负方向 $\vec{p}_i = - \nabla f(\vec{x}_i)$ 寻找，因为负梯度方向是函数值减小的最快方向。
那么不妨尝试一个例子，
\begin{eqnarray}
f(x,y) = x^2 + 10 y^2
\end{eqnarray}
显然极小值点在原点。
如果从 $(x=2, y=1)$ 出发，需要十多步才能收敛到 $(0.0062,-0.000001)$（请编程实践）。
这个方法收敛得很慢，所以前人设计了共轭梯度法。
在理论上，共轭梯度法只需要 $n$ 步，其中 $n$ 为 $\vec{x}$ 的维数；
在实践中，因为舍入误差，可能需要多几步。
相比上例中的方法，共轭梯度法是一个收敛非常快的算法。

在介绍这个算法之前，我们先介绍二次型\cite{LinearAlgebra}等概念。
\subsubsection{二次型}
若多元函数 $f(x_1, x_2, \cdots x_n)$ 有以下形式，
\begin{eqnarray}
f(x_1, x_2, \cdots, x_n) = b_{11} x^2_1 &+& b_{12} x_1 x_2 + \cdots + b_{1n} x_1 x_n \nonumber\\
																					&+& b_{22} x^2_2 + \cdots + b_{2n} x_2 x_n \nonumber\\
																					&+& \cdots + b_{nn} x^2_n
\end{eqnarray}
则称为{\bf 二次型}, 如果 $b_{ij} \in \bf R$，则 $f$ 为{\bf 实二次型}。
如果我们定义矩阵 $A$，其矩阵元为
\begin{eqnarray}
a_{ij} = \left\{
\begin{aligned}
b_{ij}, ~~~i=j,\\
\frac{1}{2} b_{ij}, ~~i \not= j,
\end{aligned}
\right.
\end{eqnarray}
就可以将 $f(\vec{x})$ 写作
\begin{eqnarray}
f(\vec{x}) = \vec{x}^\mathrm{T} A \vec{x}.
\end{eqnarray}

定理：设 $A$ 为 $n$ 阶实对称矩阵，则存在正交矩阵 $Q$，使得 $Q^{-1}AQ = Q^\mathrm{T}AQ = \Lambda$，其中
$\Lambda = \left[
\begin{array}{cccc}
\lambda_1 & & & \\
               & \lambda_2 & & \\
               &                 & \ddots &  \\
               &							&  			& \lambda_n
\end{array}
\right]$，$\lambda_1, \lambda_2, \cdots \lambda_n$ 为 $A$ 的特征值。

当 $n=1$ 时结论显然成立。
假设当 $n=k-1$ 时结论成立；设 $A$ 为 $k$ 阶实对称阵，设 $\vec{\alpha}_1$ 为 $A$ 的实特征向量，且为单位向量，$\lambda_1$ 位对应的特征值。
在线性空间中，很容易找到另外 $k-1$ 个实单位向量 $\vec{\beta}_2, \vec{\beta}_3, \cdots, \vec{\beta}_n$，使 $\vec{\alpha}_1, \vec{\beta}_2, \vec{\beta}_3, \cdots \vec{\beta}_k$ 为两两正交的单位向量组，假设 $Q_1 = (\vec{\alpha}_1, \vec{\beta}_2, \vec{\beta}_3, \cdots \vec{\beta}_k)$，则它是正交矩阵。
且
\begin{eqnarray}
Q^{-1}_1 A Q_1 &=& Q^\mathrm{T}_1 A Q_1 = \left( 
\begin{aligned}
\vec{\alpha}^\mathrm{T}_1 \\
\vec{\beta}^\mathrm{T}_2 \\
\vdots\\
\vec{\beta}^\mathrm{T}_k
\end{aligned}
\right)
A ( \vec{\alpha}_1, \vec{\beta}_2, \cdots, \vec{\beta}_k ) \nonumber\\
&=& \left(
\begin{aligned}
\lambda_1 & ~~~0 \\
0              & ~~~B_1
\end{aligned}
\right),
\end{eqnarray}
其中
\begin{eqnarray}
B_1 = \left(
\begin{aligned}
\vec{\beta}^\mathrm{T}_2 A \vec{\beta}_2 & ~~~\cdots & \vec{\beta}^\mathrm{T}_2 A \vec{\beta}_k \\
\vdots~~~~~~ & & \vdots~~~~~~\\
\vec{\beta}^\mathrm{T}_k A \vec{\beta}_2 & ~~~\cdots & \vec{\beta}^\mathrm{T}_k A \vec{\beta}_k
\end{aligned}
\right)
\end{eqnarray}
为 $k-1$ 阶实对称阵，所以存在 $k-1$ 阶正交矩阵 $P$，使得 $P^{-1} B_1 P$ 为对角阵。
令 $Q_2 = \left(
\begin{aligned}
1 & ~~~0\\
0 & ~~~P
\end{aligned}
\right)$，显然 $Q_2 A Q$ 为对角阵。
即定理得到证明。

那么，$f(\vec{x}) = \vec{x}^\mathrm{T} A \vec{x} = \vec{x}^\mathrm{T} Q^\mathrm{T} \Lambda Q \vec{x} $。
所以，如果做变换 $\vec{y} = Q \vec{x}$，即可将 $f(\vec{x})$ 写作
\begin{eqnarray}
f(\vec{x}) = \vec{y}^\mathrm{T} \Lambda \vec{y} = \lambda_1 y^2_1 + \lambda_2 y^2_2 + \cdots + \lambda_n y^2_n.
\end{eqnarray}
这种形式叫做{\bf 标准形}。
所以上述定理意味着，任何一个二次型函数，都可以化作标准形。

如果 $\lambda_1, \lambda_2, \cdots, \lambda_n$ 均大于零，则对任意 $\vec{x}$，经过相似变换 $\vec{y}=Q\vec{x}$以后，
$f(\vec{x}) = \lambda_1 y^2_1 + \lambda_2 y^2_2 + \cdots + \lambda_n y^2_n$ 一定大于等于零，当且仅当 $\vec{y}=0$ 时等于零，即当且仅当 $\vec{x}=0$时等于零（因为 $Q$ 为满秩矩阵）。
这样的实二次型叫做正定二次型，相应的矩阵 $A$ 叫做正定矩阵。


如果有一个正定对称矩阵 $A$，矢量 $\vec{p},\vec{q}$ 满足 $\vec{p}^{\mathrm{T}} A \vec{q}=0$，就说 $\vec{p},\vec{q}$ 关于 $A$ 共轭。

定理：如果 $\vec{p}_1, \vec{p}_2, \cdots, \vec{p}_n$ 中任意两个矢量均关于 $A$ 共轭，则它们必然线性无关。

%如果它们线性相关，则 $a_1 \vec{p}_1 + a_2 \vec{p}_2 + \cdots a_n \vec{p}_n =0 $ 有非零的 $a_1, a_2, \cdots, a_n$。
它们线性无关等价于下述方程有非零解 $\vec{a}$，
\begin{eqnarray}
(a_1 \vec{p}_1 + \cdots a_n \vec{p}_n)^{\mathrm{T}} A (a_1 \vec{p}_1 + \cdots + a_n \vec{p}_n) = \sum\limits^n_{i=1} a_i^2 \vec{p}_i^{\mathrm{T}} A \vec{p}_i = 0，
\end{eqnarray}
而 $A$ 为正定阵，$\vec{p} \not= 0$，即 $\vec{p}_i^{\mathrm{T}} A \vec{p}_i > 0$，所以 $a_i =0 (i=1,2,\cdots,n)$，
 $\vec{p}_1, \vec{p}_2, \cdots, \vec{p}_n$ 必然线性无关。
 这说明，这样的矢量个数 $n$ 不能大于空间维数 $N$。
 
 下面我们介绍共轭梯度法的思想。

假设 $\vec{a}$ 是 $f(\vec{x})$ 的极小值点，根据公式（\ref{multi-variable-taylor-series}），将多元函数 $f(\vec{x})$ 在 $\vec{x}=\vec{a}$ 处展开至 2 阶，得到
\begin{eqnarray}
f(\vec{x}) \approx f(\vec{a}) + \sum\limits^n_{i} f'_i (\vec{a}) (x_i-a_i) + \frac{1}{2!} \sum\limits^n_{i,j=1} f^{(2)}_{ij}(\vec{a}) (x_i - a_i)(x_j - a_j).
\end{eqnarray}
在极小值点，$f(\vec{x})$ 的 1 阶偏导数一定为零，所以
\begin{eqnarray}
f(\vec{x}) \approx f(\vec{a}) + \frac{1}{2!} \sum\limits^n_{i,j=1} f^{(2)}_{ij}(\vec{a}) (x_i - a_i)(x_j - a_j).
\end{eqnarray}
定义矩阵 $A$，其矩阵元为 $a_{ij} = \frac{1}{2} f^{(2)}_{ij} (\vec{a})$，定义 $ b = f(\vec{a})$，则有
\begin{eqnarray}
f(\vec{x}) = (\vec{x} - \vec{a})^\mathrm{T} A (\vec{x} - \vec{a}) + b,
\end{eqnarray}
注意对于一些函数形式 $f(\vec{x})$，$A$ 是有可能为零的，此处我们假设 $A$ 不是零矩阵。

由于 $\vec{a}$ 是 $f(\vec{x})$ 的极小值，所以在 $\vec{a}$ 附近的领域内，$f(\vec{x}) \geq b$。
当且仅当 $\vec{x}=\vec{a}$ 时 $f(\vec{x})=b$。
所以 A 是一个正定矩阵。
如前所述，在算法中，我们从 $\vec{x}_0$ 出发，经过 $\vec{x}_1, \vec{x}_2, \cdots$ 寻找下去，使得这些点逐步逼近极小值点 $\vec{a}$。
那么点 $\vec{x}_i$ 的梯度为
\begin{eqnarray}
\vec{g}_i \equiv \nabla f(\vec{x}_i) = 2 A (\vec{x}_i - \vec{a}).
\label{g}
\end{eqnarray}
共轭梯度方法的思路是，用 $\vec{g}_0, \cdots, \vec{g}_{i}$ 构造从 $\vec{x}_i$ 点出发寻找极小值的方向 $\vec{p}_i$（如表\ref{conjugate-gradient}），使得
\begin{eqnarray}
\forall l \not= m, ~~ \vec{g}_l \perp \vec{g}_m, ~~ \vec{p}_l^{\mathrm{T}} A \vec{p}_m=0.
\end{eqnarray}
那么，这样的 $\vec{g}_i$ 的个数不能大于空间的维数 $N$，所以理论上只需 $N$步，就可以精确找到极小值点。
当然，实际计算中存在截断误差，会多一些步数；但这个算法还是收敛得非常快，远优于前面提到的负梯度方法。
\begin{table}
\centering
\caption{\label{conjugate-gradient}共轭梯度法：给定初始值$\vec{x}_0$，沿 $\vec{p}_0$ 方向寻找函数 $f(\vec{x})$ 的极小值，记该点为 $\vec{x}_1$，求出其梯度 $\vec{g}_1$；然后沿 $\vec{p}_1$ 方向寻找极小值 $\vec{x}_2 $；如此循环下去，理论上 $ \vec{g}_n =0 $，这里 $n$ 为 $\vec{x}$ 的维数，即 $x_n$ 为极小值点。 }
\begin{tabular}{c|c|c|c}
\hline\hline
$\vec{x}_0$	&		$\vec{g}_0 = \nabla f(\vec{x}_0)$		&		$\vec{p}_0 = - \vec{g}_0	$		&		\\ \hline
$\vec{x}_1$	&		$\vec{g}_1= \nabla f(\vec{x}_1)$		&		$\vec{p}_1 = - \vec{g}_1 + \alpha_0 \vec{p}_0$	& $\alpha_0 = \frac{ \vec{g}_1^\mathrm{T} \vec{g}_1 }{ \vec{g}_0^\mathrm{T} \vec{g}_0 }$ \\
\hline
$\cdots$			&		$\cdots$				&		$\cdots$											&		$\cdots$		\\ \hline
$\vec{x}_i$	&		$\vec{g}_i=\nabla f(\vec{x}_i)$		&		$\vec{p}_i = - \vec{g}_i + \alpha_{i-1} \vec{p}_{i-1}$		& $\alpha_{i-1} = \frac{ \vec{g}_i^\mathrm{T} \vec{g}_i }{ \vec{g}_{i-1}^\mathrm{T} \vec{g}_{i-1} }$ \\
\hline
$\cdots$			&		$\cdots$				&		$\cdots$											&		$\cdots$ \\
\hline\hline
\end{tabular}
\end{table}

下面我们证明如表\ref{conjugate-gradient}所设计的思路，确实可以保证所有 $\vec{g}_0, \cdots, \vec{g}_i, \cdots$ 两两垂直，$\vec{p}_0, \cdots, \vec{p}_i, \cdots$ 两两关于 $A$ 共轭。

1. 从表\ref{conjugate-gradient}中可以看出，$\vec{p}_0, \cdots, \vec{p}_i, \cdots$ 可以用 $\vec{g}_0, \cdots, \vec{g}_i, \cdots$ 线性表示，反之亦然。这说明它们张成的线性空间完全等价。

2. 由于 $\vec{x}_1$ 是沿 $\vec{p}_0$ 方向的极小值点，所以 $\vec{g}_1 \perp \vec{p}_0$，即 $\vec{g}_1 \perp \vec{g}_0$。

3. 假设 $ \vec{g}_0, \cdots, \vec{g}_i$ 两两垂直， $\vec{p}_0, \cdots, \vec{p}_i$ 两两关于 $A$ 共轭。
如果我们可以证明 $\vec{g}_{i+1} \perp \vec{g}_0, \cdots, \vec{g}_i$，$\vec{p}_{i+1}$ 与 $\vec{p}_0, \cdots, \vec{p}_i$ 关于 $A$ 共轭，根据归纳法，我们就完成了证明。

由于 $\vec{g}_i = 2A(\vec{x}_i - \vec{a})$，所以 
\begin{eqnarray}
\vec{g}_{i+1} - \vec{g}_i = 2A(\vec{x}_{i+1} - \vec{x}_i) \propto A \vec{p}_{i} \perp \{ \vec{p}_{0}, \cdots, \vec{p}_{i-1} \}，
\end{eqnarray}
这里用到了 $\vec{p}_i$ 与 $\{ \vec{p}_0, \cdots, \vec{p}_{i-1} \}$ 关于 $A$ 共轭。
而 $\{ \vec{p}_0, \cdots, \vec{p}_{i-1} \}$ 构成的线性空间与 $\{ \vec{g}_0, \cdots, \vec{g}_{i-1} \}$ 构成的线性空间等价，所以有
\begin{eqnarray}
 \vec{g}_{i+1} - \vec{g}_i \perp \{ \vec{g}_0, \cdots, \vec{g}_{i-1} \}.
\end{eqnarray}
而 $\vec{g}_i \perp \{ \vec{g}_0, \cdots, \vec{g}_{i-1} \}$，所以有 $\vec{g}_{i+1} \perp \{ \vec{g}_0, \cdots, \vec{g}_{i-1} \}$。
需要证明的是 $\vec{g}_{i+1} \perp \vec{g}_i$。
$\vec{x}_{i+1}$ 为 $\vec{p}_i = - \vec{g}_i + \alpha_{i-1} \vec{p}_{i-1}$ 方向的极小值点，所以 
\begin{eqnarray}
\vec{g}_{i+1} \perp \vec{p}_i = -\vec{g}_i + \alpha_{i-1} \vec{p}_{i-1}.
\end{eqnarray}
$\vec{p}_{i-1}$ 由 $\{ \vec{g}_0, \cdots, \vec{g}_{i-1} \}$ 线性叠加而成， $\{ \vec{g}_0, \cdots, \vec{g}_{i-1} \}$ 均垂直于 $\vec{g}_{i+1}$，所以 $\vec{g}_{i+1} \perp \vec{g}_i$。
如此，我们证明了 $\vec{g}_{i+1} \perp \{\vec{g}_0, \cdots, \vec{g}_i\}$。

要证明 $\vec{p}_{i+1}$ 与 $\{ \vec{p}_0, \cdots \vec{p}_i \}$ 均关于 $A$ 共轭，等价于证明 
\begin{eqnarray}
\vec{p}_{i+1} = -\vec{g}_{i+1} + \alpha_{i} \vec{p}_i \perp \{ \vec{g}_1 - \vec{g}_0, \vec{g}_2 - \vec{g}_1, \cdots, \vec{g}_{i+1} - \vec{g}_i \}.
\end{eqnarray}
由于 $\vec{p}_i$ 与 $\{ \vec{p}_0, \cdots \vec{p}_{i-1} \}$ 均关于 $A$ 共轭，所以有
\begin{eqnarray}
\vec{p}_{i} \perp \{ \vec{g}_1 - \vec{g}_0, \vec{g}_2 - \vec{g}_1, \cdots, \vec{g}_{i} - \vec{g}_{i-1} \}.
\end{eqnarray}
显然 $\vec{g}_{i+1} \perp \{ \vec{g}_1 - \vec{g}_0, \vec{g}_2 - \vec{g}_1, \cdots, \vec{g}_{i} - \vec{g}_{i-1} \}$，
所以我们只需证明 
\begin{eqnarray}
\vec{p}_{i+1} = -\vec{g}_{i+1} + \alpha_{i} \vec{p}_i = (-\vec{g}_{i+1} - \alpha_i \vec{g}_i - \alpha_i \alpha_{i-1} \vec{p}_{i-1}) \perp (\vec{g}_{i+1} - \vec{g}),
\end{eqnarray}
即
$(-\vec{g}_{i+1} - \alpha_i \vec{g}_i) \perp (\vec{g}_{i+1} - \vec{g}_i)$，即
\begin{eqnarray}
(-\vec{g}_{i+1} - \alpha_i \vec{g}_i)^\mathrm{T} (\vec{g}_{i+1} - \vec{g}_i) = - \vec{g}^\mathrm{T}_{i+1} \vec{g}_{i+1} + \alpha_i \vec{g}^\mathrm{T}_i \vec{g}_i =0.
\end{eqnarray}
显然 $\alpha_i$ 的值满足要求。所以命题得到了证明。

\iffalse
%这段重写了一遍，所以注掉了。重写以后可能仍然不是特别理想。。。2017/11/14
假设从 $\vec{x}_0$ 处出发，经过数次转折，来到了 $\vec{x}_i$ 点，由于公式（\ref{g}），。
\begin{eqnarray}
\forall l = 0, \cdots, i, ~~~~ \vec{g}(\vec{x}_l) - \vec{g}(\vec{x}_{l-1}) = 2A ( \vec{x}_l - \vec{x}_{l-1} ) \propto A \vec{p}_{l-1},
\end{eqnarray}
所以
\begin{eqnarray}
\forall l=0, \cdots, i, ~~~~ \vec{g}_l - \vec{g}_{l-1} \perp \vec{p}_m, ~~~m=0, \cdots, l-1.
\label{g-g}
\end{eqnarray}
任意 $\vec{p}_m$ 只能用 $\vec{g}_s,  ~~~ s=0, \cdots, m$ 这些矢量来构造，
由于 $\vec{g}_s$ 两两垂直，所以
$\vec{p}_0, \vec{p}_1, \cdots, \vec{p}_m $
构成一个 $m+1$ 维线性空间，等价于
$\vec{g}_0, \vec{g}_1, \cdots, \vec{g}_m $ 
张成的线性空间。
公式（\ref{g-g}）意味着 $\vec{g}(\vec{x}_i) \perp \vec{g}(\vec{x}_l), ~~~~l=0,\cdots,i-1$ 是自然满足的。

所以只需构造 $\vec{p}_i$，使得它与已有的方向矢量关于 A 共轭。
假设
\begin{eqnarray}
\vec{p}_i = - \vec{g}_i + \alpha_{i-1} \vec{p}_{i-1},
\end{eqnarray}
要使它与 $\vec{p}_l, ~~~~l=0,\cdots,i-1$ 关于 A 共轭，即要求它与 $\vec{g}_1 - \vec{g}_0, \cdots, \vec{g}_i - \vec{g}_{i-1}$ 均垂直。
而 $\vec{g}_i$ 显然垂直于这些矢量，$\vec{p}_{i-1}$ 垂直于 $\vec{g}_1 - \vec{g}_0, \cdots, \vec{g}_{i-1} - \vec{g}_{i-2}$。
所以唯一的要求是 $\vec{p}_i \perp \vec{g}_i - \vec{g}_{i-1}$，即
\begin{eqnarray}
\vec{p}_i^\mathrm{T} (\vec{g}_i - \vec{g}_{i-1}) = ( - \vec{g}_i + \alpha_{i-1} \vec{p}_{i-1})^\mathrm{T} (\vec{g}_i - \vec{g}_{i-1})
= - \vec{g}_i^{\mathrm{T}} \vec{g}_i - \alpha_{i-1} ( -\vec{g}_{i-1} + \alpha_{i-2} \vec{p}_{i-2})^\mathrm{T} \vec{g}_{i-1} =0,
\end{eqnarray}
解得
\begin{eqnarray}
\alpha_{i-1} \equiv \frac{ \vec{g}_i^{\mathrm{T}} \vec{g}_i}{ \vec{g}_{i-1}^{\mathrm{T}} \vec{g}_{i-1}}.
\end{eqnarray}

这样就得到了共轭梯度法，选取参数空间中的初始位置 $\vec{x}_0$，令 $\vec{p}_0 = - \vec{g}_0$。
到达 $\vec{x}_1$后，取 $\vec{p}_1 = - \vec{g}_1 + \frac{ \vec{g}_1^\mathrm{T} \vec{g}_1 }{ \vec{g}_0^\mathrm{T} \vec{g}_0} \vec{p}_0\cdots$，
如此进行下去，对于正定二次型多元函数，理论上只需 $N$ 步即可到达极小值。
\fi

\section{第2部分作业}
１．利用Richardson算法，计算$e^x$在$x=10$处的二阶导数，计算$D^0_3$，并计算误差。

*2. 牛顿也写了一个多项式插值公式，他定义
\begin{eqnarray}
P_n (x) &=& a_0 + a_1 (x-x_0) + a_2(x-x_0)(x-x_1) \nonumber\\
        && + \cdots + a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1}),
\end{eqnarray}
它的前1项通过第一个数据点，前2项通过前两个数据点，...，前n+1项在一起
通过前n+1个数据点 $(x_0,y_0), (x_1,y_1), \cdots, (x_n,y_n)$。

1)试证明
\begin{equation}
a_n = \frac{ y_n - P_{n-1} (x_n) }{ (x_n-x_0)(x_n-x_1)\cdots(x_n-x_{n-1}) },
\end{equation}

2)试用递归函数编写一个程序，自动读入数据，用牛顿插值公式进行拟合，并给出指定x值处的插值$P_n(x)$。

3．粒子从$t=0s$到$t=10s$，速度为$v(t)=10t-t^2 (m/s)$，请用梯形法，积分得到粒子走过的路程，并与解析解相比较，误差为多大？

*4.　编写10阶高斯-勒让德积分，计算$\int^1_0 e^x dx$，并与理论值比较，得出误差。

5.　将$[0,\pi]$分割为10000份，用中值法、梯形法、辛普生方法分别积分$\int^\pi_0 \sin(x) dx$，比较误差，并进行简单分析。

*6.　试讨论：怎样构造自适应算法+梯形法进行积分？

7.试选取适当的区间，找到$f(x)=x(x-1)(x-2)(x-3)$的两个极小值。

*8.寻找$f(x)=x(x-1)(x-2)(x-3)$的一个极大值。

9.试编写辛普生积分的openmp并行程序，并自选一个任务，记录并行v.s.单线程完成任务所需时间。

10.任选一种方法，找到$f(x) = x^2 -4x -5$的两个零点，要求精度在$10^{-9}$以内，并给出迭代次数。

*11.试用Stewenson方法，找到$f(x) = x^2 -4x -5$的两个零点，要求精度在$10^{-9}$以内，并给出迭代次数。